      SUBROUTINE E01AEY(WITHQ0,M,X,XMIN,XMAX,Y,IP,IMAX,N,NP1,ITMIN,
     *                  ITMAX,A,LA,RES,PMAX,PINDEX,NIT,ATRIAL,PTRIAL,
     *                  FTAU,C,D,W,ADIF,DA,RNM,RTRLNM,LOCX,LOCY,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE REFH
C
C     CREATED 17 07 79.  UPDATED 14 05 80.  RELEASE 00/47
C
C     AUTHORS ... GERALD T. ANTHONY, MAURICE G. COX
C                 J. GEOFFREY HAYES AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEY.  A ROUTINE TO APPROXIMATE AND THEN REFINE AN
C     INTERPOLATING POLYNOMIAL  Q(X)  OR A ZEROIZING
C     POLYNOMIAL  Q0(X)  IN ITS CHEBYSHEV REPRESENTATION
C
C     INPUT PARAMETERS
C        WITHQ0   TRUE IF ZEROIZING POLYNOMIAL, ELSE FALSE
C        M        THE NUMBER OF DISTINCT DATA POINTS.
C        X        ARRAY CONTAINING THE DISTINCT X-VALUES,
C                    NORMALIZED IF NECESSARY TO (-1, 1).
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C      * Y        ARRAY CONTAINING VALUES AND DERIVATIVES OF
C                    THE DEPENDENT VARIABLE.
C        IP       ARRAY SPECIFYING THE HIGHEST ORDER OF
C                    DERIVATIVE AT EACH X-VALUE.
C        IMAX     ONE MORE THAN THE LARGEST ELEMENT OF THE
C                    ARRAY IP.
C        N        NUMBER OF INTERPOLATING CONDITIONS.
C                    N = M + IP(1) + IP(2) + ... + IP(M).
C        NP1      VALUE OF  N + 1
C        ITMIN,
C        ITMAX    THE LOWER AND UPPER LIMITS ON THE ITERATIVE
C                    PROCESS.
C
C     OUTPUT (AND ASSOCIATED DIMENSION) PARAMETERS
C        A        CHEBYSHEV COEFFICIENTS OF POLYNOMIAL
C        LA       DIMENSION OF  A.
C                    .GE. N  IF INTERPOLATING POLYNOMIAL
C                    .GE. N + 1  IF ZEROIZING POLYNOMIAL
C        RES      RESIDUALS OF POLYNOMIAL
C        PMAX     LARGEST PERFORMANCE INDEX
C        PINDEX   PERFORMANCE INDICES
C        NIT      NUMBER OF ITERATIONS TAKEN
C
C     WORKSPACE PARAMETERS
C        ATRIAL   TRIAL VALUES OF THE CHEBYSHEV COEFFICIENTS
C        PTRIAL   PERFORMANCE INDICES CORRESPONDING TO  ATRIAL
C      * FTAU     SCALED VALUES OF  Y.  IF  Y(I)  IS THE
C                    VALUE OF AN  R-TH  DERIVATIVE, THEN
C                    ((XMAX - XMIN)/2)**R/(FACTORIAL R)
C                    TIMES  Y(I)  IS THE VALUE OF  FTAU(I)
C      * C        COEFFICIENTS IN NEWTON FORM OF POLYNOMIAL
C      * D        INTERMEDIATE DIVIDED DIFFERENCE VALUES
C      **W        VALUES OF CORRECTION POLYNOMIAL AT
C                    CHEBYSHEV EXTREMA
C        ADIF     CHEBYSHEV COEFFICIENTS OF A DERIVATIVE OF
C                    AN APPROXIMATION TO THE POLYNOMIAL
C        DA       CHEBYSHEV COEFFICIENTS OF A
C                    CORRECTION POLYNOMIAL
C        RNM      RESIDUAL NORMS CORRESPONDING TO  A
C        RTRLNM   RESIDUAL NORMS CORRESPONDING TO  ATRIAL
C      * LOCX     POINTERS TO X-VALUES IN CONSTRUCTING
C                    NEWTON FORM OF POLYNOMIAL
C      * LOCY     POINTERS TO Y-VALUES CORRESPONDING TO X-VALUES
C
C     FAILURE INDICATOR PARAMETER
C        IFAIL    FAILURE INDICATOR
C                    0 - SUCCESSFUL TERMINATION
C                    1 - ITERATION LIMIT EXCEEDED
C                    2 - ITERATION DIVERGENT
C
C           NOTES.  (1) THE ELEMENTS OF THE ARRAYS MARKED  *  ARE
C                       NOT ACCESSED IF  WITHQ0  IS  TRUE.
C                   (2) THE ELEMENTS OF THE ARRAY MARKED  **  IS
C                       NOT ACCESSED IF  WITHQ0  IS  FALSE.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  PMAX, XMAX, XMIN
      INTEGER           IFAIL, IMAX, ITMAX, ITMIN, LA, M, N, NIT, NP1
      LOGICAL           WITHQ0
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), ADIF(LA), ATRIAL(LA), C(N), D(N), DA(LA),
     *                  FTAU(N), PINDEX(IMAX), PTRIAL(IMAX), RES(N),
     *                  RNM(IMAX), RTRLNM(IMAX), W(NP1), X(M), Y(N)
      INTEGER           IP(M), LOCX(M), LOCY(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  AMAX, ATRNRM, DANRM, HALF, ONE, PMXTRL, SCALE,
     *                  SXTEEN, ZERO
      INTEGER           I, IERROR, IT, ITEMP, ITMXP1, ITP1, L, NFREF,
     *                  NPILT1, NTERMS
      LOGICAL           IMPROV, WITHPI, ZERODA
C     .. External Subroutines ..
      EXTERNAL          E01AEU, E01AEV, E01AEZ
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, LOG
C     .. Data statements ..
      DATA              ZERO, HALF, ONE, SXTEEN/0.0D+0, 0.5D+0, 1.0D+0,
     *                  16.0D+0/
C     .. Executable Statements ..
      IERROR = 0
C
C     NUMBER OF TERMS IN POLYNOMIAL
C
      NTERMS = N
      IF (WITHQ0) NTERMS = NTERMS + 1
C
C     NUMBER OF PERFORMANCE INDICES LESS THAN ONE
C
      NPILT1 = 0
C
C     INDICATE THAT PERFORMANCE INDICES ARE TO BE PRODUCED
C     ONLY IF AN IMPROVEMENT IS OBTAINED
C
      WITHPI = .FALSE.
C
C     INDICATE THAT THE FINE REFINEMENT STAGE HAS NOT YET STARTED
C
      NFREF = -2
C
C     SET RESIDUALS INITIALLY EQUAL TO SPECIFIED Y-VALUES
C     (IF  Q(X)  REQUIRED) OR ZERO (IF  Q0(X)  REQUIRED)
C
      DO 20 I = 1, N
         IF ( .NOT. WITHQ0) RES(I) = Y(I)
         IF (WITHQ0) RES(I) = ZERO
   20 CONTINUE
C
C     INITIALIZE TRIAL CHEBYSHEV COEFFICIENTS
C
      DO 40 I = 1, NTERMS
         ATRIAL(I) = ZERO
   40 CONTINUE
C
C     COMMENCE ITERATIVE REFINEMENT
C
      ITMXP1 = ITMAX + 1
      DO 360 ITP1 = 1, ITMXP1
C
C        IT  IS THE ACTUAL ITERATION NUMBER,  IT = 0
C        CORRESPONDING TO THE FIRST ESTIMATE OF THE POLYNOMIAL
C
         IT = ITP1 - 1
C
C        DETERMINE CHEBYSHEV COEFFICIENTS  DA  OF POLYNOMIAL
C        APPROXIMATELY SATISFYING THE CONDITIONS IN  RES
C
         IF (WITHQ0 .AND. IT.EQ.0) CALL E01AEU(M,X,IP,NP1,DA,W)
         IF ( .NOT. WITHQ0 .OR. (WITHQ0 .AND. IT.GT.0))
     *       CALL E01AEZ(M,XMIN,XMAX,X,RES,IP,N,DA,LOCX,LOCY,FTAU,D,C)
C
C        SKIP TEST FOR DIVERGENCE IF ON ZERO-TH ITERATION
C
         IF (IT.EQ.0) GO TO 100
C
C        DETERMINE THE NORMS OF  DA  AND (THE PREVIOUS)  ATRIAL
C
         DANRM = HALF*ABS(DA(1))
         ATRNRM = HALF*ABS(ATRIAL(1))
         IF (N.EQ.1) GO TO 80
         DO 60 I = 2, N
            DANRM = DANRM + ABS(DA(I))
            ATRNRM = ATRNRM + ABS(ATRIAL(I))
   60    CONTINUE
   80    IF (WITHQ0) ATRNRM = ATRNRM + ABS(ATRIAL(NP1))
C
C        ASSUME DIVERGENCE IF THE NORM OF  DA  IS NOT
C        LESS THAN THAT OF  ATRIAL  ...
C
         IF (DANRM.GE.ATRNRM) IERROR = 2
         IF (DANRM.GE.ATRNRM) GO TO 380
C
C        ... OTHERWISE DETERMINE NEW TRIAL APPROXIMATION
C
  100    ZERODA = .TRUE.
         DO 120 I = 1, N
            ATRIAL(I) = ATRIAL(I) + DA(I)
            IF (DA(I).NE.ZERO) ZERODA = .FALSE.
  120    CONTINUE
         IF ( .NOT. (WITHQ0 .AND. IT.EQ.0)) GO TO 140
         ATRIAL(NP1) = DA(NP1)
         IF (DA(NP1).NE.ZERO) ZERODA = .FALSE.
C
C        DETERMINE RESIDUALS, PERFORMANCE INDICES AND
C        LARGEST PERFORMANCE INDEX CORRESPONDING TO
C        TRIAL COEFFICIENTS  ATRIAL
C
  140    CALL E01AEV(WITHQ0,WITHPI,M,XMIN,XMAX,X,N,Y,IP,IMAX,ATRIAL,LA,
     *               IT,RNM,RTRLNM,IMPROV,ADIF,RES,PMXTRL,PTRIAL)
C
C        SET DUMMY, NON-ZERO, VALUE OF  PMXTRL  IF NO
C        IMPROVEMENT, OTHERWISE IT IS UNDEFINED
C
         IF ( .NOT. IMPROV) PMXTRL = SXTEEN
C
C        IF ON FIRST ITERATION, OR IF THE LARGEST PERFORMANCE
C        INDEX IS ZERO, OR IF ALL COMPONENTS OF  DA  ARE ZERO,
C        TAKE THE TRIAL SET OF COEFFICIENTS AND PERFORMANCE
C        INDICES AS THE BEST (SO FAR)
C
         IF ( .NOT. (IT.EQ.0 .OR. PMXTRL.EQ.ZERO .OR. ZERODA))
     *       GO TO 200
         DO 160 I = 1, NTERMS
            A(I) = ATRIAL(I)
  160    CONTINUE
         DO 180 L = 1, IMAX
            RNM(L) = RTRLNM(L)
            PINDEX(L) = PTRIAL(L)
  180    CONTINUE
         PMAX = PMXTRL
C
C        FINISH IF LARGEST PERFORMANCE INDEX IS ZERO
C        OR IF ALL COMPONENTS OF  DA  ARE ZERO
C        (I.E. NO FURTHER IMPROVEMENT IS POSSIBLE)
C
  200    IF (PMXTRL.EQ.ZERO .OR. ZERODA) GO TO 380
C
C        INDICATE WHETHER THE FINE REFINEMENT STAGE HAS COMMENCED
C        (I.E. FOR THE FIRST TIME ALL PERFORMANCE INDICES ARE
C        LESS THAN ONE)
C
         IF (NFREF.EQ.-2 .AND. PMXTRL.LT.ONE) NFREF = -1
C
C        BRANCH ACCORDING TO WHETHER THE PROCESS IS IN THE
C        FINE REFINEMENT STAGE  (NFREF .GE. 0)  OR NOT
C        (NFREF .EQ. -1)
C
         IF (NFREF.GE.-1) GO TO 280
C
C        THE PROCESS IS IN THE COURSE REFINEMENT PHASE.
C        UPDATE THE COEFFICIENTS AND THE CORRESPONDING
C        NORMS AND PERFORMANCE INDICES IF
C           (I)  THERE HAS BEEN AN IMPROVEMENT IN (AT
C                LEAST) ONE OF THE RESIDUAL NORMS, AND
C           (II) THE NUMBER OF PERFORMANCE INDICES
C                THAT ARE LESS THAN ONE HAS NOT
C                INCREASED COMPARED WITH THOSE OF
C                THE BEST POLYNOMIAL SO FAR.
C
         IF ( .NOT. IMPROV) GO TO 360
         ITEMP = 0
         DO 220 L = 1, IMAX
            IF (PTRIAL(L).LT.ONE) ITEMP = ITEMP + 1
  220    CONTINUE
         IF (ITEMP.LT.NPILT1) GO TO 360
         NPILT1 = ITEMP
         DO 240 I = 1, NTERMS
            A(I) = ATRIAL(I)
  240    CONTINUE
         DO 260 L = 1, IMAX
            RNM(L) = RTRLNM(L)
            PINDEX(L) = PTRIAL(L)
  260    CONTINUE
         PMAX = PMXTRL
         GO TO 360
C
C        THE PROCESS IS IN THE FINE REFINEMENT PHASE.
C        UPDATE THE COEFFICIENTS AND THE CORRESPONDING
C        NORMS AND PERFORMANCE INDICES IF
C           (I)  THERE HAS BEEN AN IMPROVEMENT IN (AT
C                LEAST) ONE OF THE RESIDUAL NORMS, AND
C           (II) THE LARGEST PERFORMANCE INDEX IS LESS
C                THAN THE LARGEST OF THAT OF THE BEST
C                POLYNOMIAL SO FAR.
C        INCREMENT THE NUMBER OF FINE REFINEMENTS (THE NUMBER
C        OF REFINEMENTS SINCE THE FIRST OCCASION WHEN ALL
C        PERFORMANCE INDICES WERE LESS THAN UNITY), EXITING
C        IF AS MANY AS  ITMIN  FINE REFINEMENTS HAVE BEEN
C        PERFORMED
C
  280    IF ( .NOT. IMPROV) GO TO 340
         IF (PMXTRL.GE.PMAX) GO TO 340
         DO 300 I = 1, NTERMS
            A(I) = ATRIAL(I)
  300    CONTINUE
         DO 320 L = 1, IMAX
            RNM(L) = RTRLNM(L)
            PINDEX(L) = PTRIAL(L)
  320    CONTINUE
         PMAX = PMXTRL
  340    NFREF = NFREF + 1
         IF (NFREF.GE.ITMIN) GO TO 380
  360 CONTINUE
C
C     THE PROCESS HAS NOT SUCCEEDED IN REDUCING ALL THE
C     PERFORMANCE INDICES TO LESS THAN UNITY
C
      IERROR = 1
      IT = ITMAX
C
C     NUMBER OF ITERATIONS ACTUALLY TAKEN
C
  380 NIT = IT
      IF ( .NOT. WITHQ0) GO TO 440
C
C     IN THE CASE OF  Q0(X),  SCALE ITS COEFFICIENTS BY
C     AN INTEGRAL POWER OF  16  SUCH THAT THE LARGEST
C     COEFFICIENT IS OF ORDER UNITY
C
      AMAX = ZERO
      DO 400 I = 1, NP1
         IF (ABS(A(I)).GT.AMAX) AMAX = ABS(A(I))
  400 CONTINUE
      IF (AMAX.EQ.ZERO) GO TO 440
      I = LOG(AMAX)/LOG(SXTEEN)
      SCALE = SXTEEN**(-I)
      DO 420 I = 1, NP1
         A(I) = SCALE*A(I)
  420 CONTINUE
C
C     RETURN RESIDUALS CORRESPONDING TO SELECTED COEFFICIENTS
C
  440 WITHPI = .TRUE.
      CALL E01AEV(WITHQ0,WITHPI,M,XMIN,XMAX,X,N,Y,IP,IMAX,A,LA,IT,RNM,
     *            RTRLNM,IMPROV,ADIF,RES,PMAX,PINDEX)
      IFAIL = IERROR
      RETURN
C
C     END E01AEY
C
      END
